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Abstract. Algebraic multigrid is an iterative method that is often optimal for solving the matrix equations that arise 
in a wide variety of applications, including discretized partial differential equations. It automatically constructs a sequence 
of increasingly smaller matrix problems that enable efficient resolution of all scales present in the solution. One of the main 
components of the method is an adequate choice of coarse grids. The current coarsening methodology is based on measuring how 
a so-called algebraically smooth error value at one point depends on the error values at its neighbors. Such a concept of strength 
of connection is well understood for operators whose principal part is an M-matrix; however, the strength concept for more 
general matrices is not yet clearly understood, and this lack of knowledge limits the scope of AMG applicability. The purpose 
of this paper is to motivate a general definition of strength of connection, based on the notion of algebraic distances, discuss its 
implementation, and present the results of initial numerical experiments. The algebraic distance measure, we propose, uses as 
its main tool a least squares functional, which is also applied to define interpolation. 
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1. Introduction. Multigrid methods for solving general systems of linear equations Ax — b are all 
based on the smoothing property of relaxation. For any < r < 1, an error vector e is called r-smooth if 
all its residuals are smaller than r||e||. The basic observation Bra86 is that the convergence of a proper 
relaxation process slows down only when the current error is r-smooth with r <C 1, the smaller the r the 
slower the convergence. This observation and the assumption that when relaxation slows down, the error 
vector e can be approximated in a much lower-dimensional subspace, are the main ideas behind the multigrid 
methodology. Very efficient geometric multigrid solvers have been developed for the case that the lower- 
dimensional subspace corresponds to functions on a well-structured grid (the coarse level), on which the 
smooth errors can be approximated by easy-to-derive equations, based for example on discretizing the same 
continuum operator that has given rise to the fine-level equations Ax — b. The coarse-level equations are 
solved recursively using a similar combination of relaxation sweeps and still-coarser-level approximations to 
the resulting smooth errors. 

To deal with more general problems, e.g., ones where the fine- level system may not be defined on a 
well-structured grid nor perhaps arise from any continuum problem, algebraic multigrid (AMG) methods 
develop techniques for deriving the set of coarse-level variables and the coarse-level equations, based directly 
on the given matrix A. The basic approach, developed in BMR83 ( BMR84 ( RS87 and called today classical 
AMG or RS-AMG, involves the following two steps. 

(1) Choosing the coarse-level variables as a subset C of the set O of fine variables, such that each variable 
in fi is strongly connected to variables in C. 

(2) Approximating the fine-level residual problem Ae — r by the coarse-level equations A c e c — r c using 
the Galerkin prescription A c = RAP and r c = Rr, yielding an approximation Pe c to e. 

The interpolation matrix P, and the restriction matrix R are both defined directly in terms of the elements 
of the matrix A, and this construction relies heavily on the notion of strong connections developed for 
M-matrices. 

In the past two decades, many variations of the classical AMG algorithm have been introduced, including 
modifications to the coarse-grid selection algorithms CFHJ98 DYH06 HY02 , the definition of interpola- 
tion 



BFM+05 , or both BL07 BCF+00 CFH+03 VMB96 



variations within the classical AMG framework appears in MeuOl 



A summary and comparison of many of these 
Such variations arise because the ap- 



plicability of the original algorithm and its variants is limited by the Ai-matrix heuristics upon which they 
are based. In particular, these variants of the classical AMG approach yield efficient solvers for problems 
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where the (properly scaled) matrix A has a dominant diagonal and, with small possible exceptions, all its 
off-diagonal elements have the same sign. Even then, the produced interpolation can have limited accuracy, 
insufficient for full multigrid efficiency. An example where performance of the classical AMG method can 
deteriorate, and the one we consider here, is given by scalar anisotropic diffusion problems. 

Recent advances in the development of algebraic multigrid have focused on alternative notions of strength 



of connection that do not rely on specific properties of the system matrix, A, in their construction. In Cho03 



for example, samples of algebraically smooth error are used to determine the directions in which this error 
varies slowly and, thus, interpolation can be very effective. For many problems, however, such a direction 
may not exist, and the proposed technique does not apply. An alternative approach was developed by 
Broker Bro03] , in which the relative sizes of the entries of the sparse approximate inverse preconditioning 



matrix |GH97 are used to indicate strength of connection. In some respects, this approach is a special case 



of the approach proposed in BBM+06 , which also measures strength using columns of the inverse; however, 



the measure used in |Bro03 is based on the Li sizes of the approximate inverse, not on the energy norm 
as in |BBM + 06 . Closely related to this work is the approach of compatible relaxation Bra00,BF10 BL07 



FV04 Liv04 which uses a modified relaxation scheme to expose the character of the slow-to-converge, i.e., 



algebraically smooth, error. Coarse-grid points are then selected where this error is the largest, thus avoiding 
explicit use of a measure of strength of connection in choosing the coarse variables. 

In the present article, we focus on developing an alternative coarsening algorithm to select coarse variable 
and interpolatory sets for anisotropic systems where even these newly developed variants of AMG can often 
be difficult to apply. Our approach is based on a general notion of strength of connection derived from 
the minimum of the local least squares problem used in defining interpolation [BraOl BBKL11 . The basic 



approach finds interpolation to fit a set of test vectors (representatives of algebraically smooth error) in a 
least squares sense. We use the LS functional resulting from a simplified direct form of interpolation to 
construct an algebraic distance measure of strength of connection which is used to derive a strength graph. 

Chap. 8] to coarsen the unknowns at 



BHMOO 



The strength graph is then passed to a coloring algorithm 
each stage of a compatible relaxation coarsening algorithm. The quality of interpolatory sets are measured 
a posteriori, by first constructing LS interpolation to these sets and then monitoring the corresponding 
values of the LS functional for the computed minimizers BraOl BBKL11 . A main new feature of our 



coarsening algorithm is that it allows for aggressive and problem-oriented coarsening as well as long-range 
interpolation in a straightforward way. This is accomplished by deriving long-range measures of strength 
between neighboring pairs (or sets) of fine and coarse points. 

Several algorithms for aggressive coarsening and long-range interpolation have been developed in the 
past Stu01||DYH06 , yet these techniques rely on the same M-matrix based strength of connection measure 
used in the original classical AMG coarsening algorithm BMR83 BMR84, RS87 . Aggressive coarsening is 
then achieved by applying the classical AMG coloring algorithm multiple times before constructing inter- 
polation. While such approaches have been shown to be effective in reducing the large complexities often 
resulting from application of the classical AMG scheme to finite-difference and finite-clement discretizations 
of elliptic problems, they are again limited to M-matrix type problems. In contrast, our proposed long-range 
measure aims to address the issue of strength of connection in a more general context the goal being to 
determine explicitly those degrees of freedom from which high quality least squares interpolation for some 
given set of test vectors can be constructed. 

The remainder of the paper is organized as follows. First, in Section [2j we give an overview of the 
classical algebraic multigrid coarsening approach. Section [3] contains an introduction to the bootstrap alge- 
braic multigrid components, with an emphasis on the least squares interpolation approach and compatible 
relaxation coarsening algorithm. Then, in Section |4j a general definition of strength of connection and the 
notion of algebraic distances, as well as its connection to compatible relaxation, are discussed. The diffusion 
equation with anisotropic coefficients and its descritization are introduced in Section [5] as are the results of 
numerical experiments of our approach applied to this system. 

2. Classical AMG Coarsening. The coarse-grid selection performed in the classical AMG algorithm 



BMR84 , RS87 is made based on properties of slow-to-converge errors discovered by applying pointwise 
relaxation to discrete systems with diagonally dominant M-matrices. These errors are used to identify 
the important (strong) connections within the linear system. The coarse-grid points are then selected using 
maximal independent subset heuristics to ensure a significant reduction in grid size, still maintaining accurate 
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approximation properties. The outline of the classical AMG coarsening procedure is given next. 

Consider the Gauss-Seidel iteration, with error propagation operator I — L~ 1 A 1 where L is the lower- 
triangular part of A. Fixing this to be the fine-scale relaxation used in a multigrid method for the symmetric 
and diagonally dominant M-matrix, A, the goal of coarse-grid correction is to effectively reduced the error 



components that are not significantly reduced by Gauss-Seidel relaxation. It is shown in BMR84| that such 
an error, e, must satisfy 



Thus, if etjj is large, relative to max \a,ik\ or max \cijk\t then it must be true that rs ej. 

This observation leads to the definitions of strong dependence which has significant influence on the 
selection of AMG coarse-grid points. For a given degree of freedom, i, the set of points that i strongly 
depends upon is defined as 



Si = < j : -dij > 0m&x{-a ik } > , (2.1) 
L k ^ J 

for some suitable choice of 6, < < 1. Similarly, the set of points that are strongly influenced by i is 
denoted Sf: j £ Sf implies that i £ Sj. 

Once strong connections are determined, a coarse grid is chosen so that all strongly connected neighbors 
of any fine-grid point (that is not itself a coarse-grid point) are available for direct or a path-length two 
indirect interpolation. That is, if two fine-grid points are strongly connected to one another, they must have 
a common coarse-grid neighbor, so that this strong connection may be accounted for indirectly in interpo- 
lation. This condition is usually expressed as the first AMG coarsening heuristic. 

HI: For each fine-grid point i, each point j £ Si must either be a coarse-grid neighbor of i or strongly 
depend on at least one coarse-grid neighbor of i. 

If used alone, this rule tends to create large coarse grids, as it is most easily satisfied by adding points to 
the coarse grid. In practice, a second heuristic is used to limit the size of the coarse grid. 

H2: The set of coarse-grid points, C, should be a maximal subset of the original set of fine-grid points 
such that no coarse-grid point strongly depends on any other coarse-grid point. 

Because these two goals may be contradictory, typical AMG implementations rely on enforcing HI and use 
H2 as a guide. This is done by selecting coarse-grid points using a two-pass algorithm that picks a set 
satisfying H2, then checks for any points where HI is violated, adding new coarse-grid points to compensate 



if this occurs. The first stage is often implemented as a coloring algorithm BHMOO Chap. 8], where coarse 
points are selected based on their number of strongly connected neighbors. Initially, all points are weighted 
with the number of points that strongly depend upon them (that is, the size of Sf). The point with the 
largest weight is then selected to be a coarse-grid point. Since each j £ Sf is now strongly connected to 
a coarse-grid point, all such j are made fine-grid points so that H2 is not violated. All strongly connected 
neighbors of these points (that is, k £ Sj for any j £ Sf) are then made more attractive as coarse-grid 
points, since they reflect unresolved strong connections of fine-grid points. Thus, the weights of all such k 
are incremented for each j £ Sk that was made a fine-grid point. The algorithm then repeats, selecting the 
new largest-weighted point as a coarse-grid point. 

In this way, an initial coarse grid is chosen to give a maximal independent set over all strong connections. 
Then, if necessary, the minimal number of points needed to ensure that property HI is fulfilled are added to 
the coarse grid. Once the selection of coarse grid has been completed, an interpolation operator is defined 
so that it is accurate for errors that are slowly varying along strong connections. We do not provide specific 



details of AMG interpolation here and, instead, again refer the reader to BHMOO Chap. 8]. We note, 
however, that the classical approach for constructing P also relies on the above M-matrix heuristic that the 
direction in which the smooth error is locally constant can be identified by the entries in the system matrix, 
thus, limiting its range of applicability. 
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3. Bootstrap AMG. The bootstrap AMG (BAMG) algorithm, introduced in [Bra01[ §17.2], was 
developed to extend AMG to general (non M-matrix) problems. The BAMG approach combines the following 
two general devices to inexpensively construct high-quality interpolation. 

(A) The interpolation is derived to provide the best least squares fit to a set of r-smooth test vectors 
(TVs) obtained by a process described below. 

Denote by Cj the set of coarse variables used in interpolating to i £ F. It follows from the satisfaction of 
the compatible relaxation criterion (see the next section) that, with a proper choice of Ci, there exists an 
interpolation that will commit only 0(r) errors in reproducing Xi for any vector x which is r-smooth. The 
size \Ci\ of this set should in principle increase as r decreases (and smaller r means overall better multigrid 
convergence), but in practice a pre-chosen and sufficiently small interpolation caliber c := maxjgjr |C,-|, often 
yields small enough errors. The set Ci can often be adequately chosen by natural considerations, such 
as the set of geometrical neighbors with i in their convex hull. If the chosen set is inadequate, the least 
squares procedure will show bad fitness (interpolation errors large compared with r), and the set can then 
be improved. The least squares procedure can also detect variables in d that can be discarded without 
significant accuracy loss. Thus, this approach allows creating interpolation with whatever needed accuracy 
which is as sparse as possible. 

(B) Generally, the test vectors are constructed in a bootstrap manner, in which several tentative AMG 
levels are generated by interpolation fitted to only moderately smooth TVs; this tentative (multi- 
level) structure is then used to produce improved (much smoother) TVs, yielding a more accurate 
interpolation operator. The process continues if needed until fully efficient AMG levels have been 
generated. 

The first test vectors are each produced by relaxing the homogeneous system Ax = with a different 
starting vector. This quickly leads to a r-smooth test set with r< 1. (A mixture of random vectors and/or 
diverse geometrically smooth vectors can generally be used as initial approximations. In the case of discretized 
isotropic PDEs, if geometrically smooth vectors that satisfy the homogeneous boundary conditions are used, 
relaxation may not be needed at all. In many other cases relaxation can be confined to the neighborhood 
of boundaries and discontinuities.) For the anisotropic problem considered, we use a constant vector and 
a set of, initially positive random, relaxed test vectors to define an algebraic distance notion of strength of 
connection and the least squares interpolation operator. 

3.1. Least squares interpolation. The basic idea of the least squares interpolation approach is to 
approximate a set of test vectors, V = {v^\ . . . ,t/ fc )} C K n , minimizing the interpolation error for these 
vectors in a least squares sense. In the context considered here, namely, applying the least squares process to 
construct a classical AMG form of interpolation, each row of P, denoted by pi, is defined as the minimizer of a 
local least squares functional. Given a splitting of variables F = Q\C for each i £ F find pi £ R" c , n c = |C|, 
such that 

where C,; C C. The weights u K > can be chosen to reflect the global algebraic smoothness of the test 
vectors. We give our specific choice in the numerical experiments section. 

Remark 3.1. We note that if the restricted vectors v c , . . . , v c ; form a basis for the local linear space 
W H , n, = \C{\, then the solutions to the local least squares minimization problems are unique. This in turn 
suggests setting a lower bound on the number of vectors, k, used in the least squares fit 

k>c. 

Further, in practice we have observed that the accuracy of the least squares interpolation operator and, hence, 
the performance of the resulting solver generally improves with increasing k \BBKLiq , up to some value 
proportional to caliber c. 

Remark 3.2. In fKahOS^ , it was observed that the implicit application of a local Jacobi relaxation to the 
TVs used in the least squares definition of interpolation is equivalent to an operator induced form of least 
squares interpolation constructed using an Element-free AMG type approach (see \Vas09l). This equivalence 
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of the least squares approach with an additional local relaxation step was also formulated and discussed in a 
slightly different scope in \MMPRlOj , where it was defined as a residual-based least squares fit 

csu, > - \ ' <" ' • ~ E &)M K) ) * min ' (3-2) 

which in turn was shown to be consistent with a classical AMG operator induced form of least squares 
interpolation. 

3.2. Compatible relaxation. A general criterion for choosing an adequate set of coarse variables 
is the fast convergence of compatible relaxation (CR), as introduced in Bra00| (being a special case of 




choosing coarse variables for much more general types of systems BralO , introduced first for problems in 
statistical mechanics [BR]). In fact, an improved version called habituated compatible relaxation, introduced 
in 



Liv04 , yields an accurate prediction for the convergence rate that can be achieved by an AMG solver that 
employs the given relaxation scheme and the proposed set of coarse variables. Analysis of general compatible 
relaxation approaches and works on the development of compatible relaxation coarsening algorithms are 



found in FV04|[BL07 BF10 . We note that the ability of habituated compatible relaxation to quantitatively 



predict performance of the solver is very useful in designing and debugging the actual solver. 

Given the coarse-and-fine level splitting, C and F, a simple form of compatible relaxation is given by 
_F-relaxation for the homogeneous system — relaxation applied only to the set of F variables. Given the 
partitioning of into F and C, we have 

U=M, B=( B Jf '!/■) . andiW < ^ 



u c J ' \B cf B CC J ' \M cf M, 

assuming the equations are permuted such that the unknowns in F come before those in C . The ^-relaxation 
form of compatible relaxation is then defined by 

= ( J - M j} B ff) u ) = (3-3) 
If M is symmetric, the asymptotic convergence rate of compatible relaxation, 

Pf = P( E f), 

where p denotes the spectral radius, provides a measure of the quality of the coarse space, that is, a measure 
of the ability of the set of coarse variables to represent error not eliminated by the given fine-level relaxation. 
This measure can be approximated using i^-relaxation for the homogeneous system with a random initial 
guess Uj. Since ]im v -^. 00 ||£^|| 1 < /,y = p(Ef) for any norm || • ||, the measure 

g f = (M\\/H\\) 1/v (3-4) 

estimates pf. 

The set of coarse variables is then constructed using a multistage coarsening algorithm, where a single 
stage consists of: (1) running several iterations of compatible relaxation (based on the current set F) and 
(2) if it is slow to converge, adding an independent set of the fine-level variables not well treated by CR to 
C. Steps (1) and (2) are applied repeatedly, giving rise to a sequence of coarse variable sets: 

= C C d C ... c C m , 

until convergence of compatible relaxation below a prescribed tolerance is reached, yielding an accepted 
coarse set C := C m . 

We note that all current CR algorithms have intentionally avoided the explicit use of a strength of 
connection measure in constructing the coarse variable sets, instead relying on the error produced by CR to 
form candidate sets of potential C-points. Our aim here is to develop a more general notion of strength of 
connections based on algebraic distances and to explore its use in a compatible relaxation coarsening scheme 
and defining AMG interpolation. 



() 



A. BRANDT, J. BRANNICK, K. KAHL, AND I. LIVSHITS 



4. Selecting coarse variables and interpolation via algebraic distances. While the classical 
definition of strong dependence is appropriate for the case of diagonally dominant M-matrices for which 
it was intended, it is frequently seen to break down when applied to problems involving other classes of 
matrices. The near null space of a diagonally dominant M-matrix is typically a locally slowly varying (or 
locally constant) vector, and, in such cases the AMG strength heuristic relies on this being reflected in the 
coefficients of the system matrix, A. If either the near null space cannot be accurately characterized as 
locally constant or this is not reflected in the matrix coefficients, then AMG performance typically suffers. 

To derive a more general measure of strength of connection, we consider qualifying strength among 
variables based on a variables ability to interpolate r-smooth error to its neighbors. We thus characterize 
a given variable (set of variables) as strongly connected to a neighboring grid point as one(s) for which the 
value of the LS functional is small when building interpolation from this point (set of points), for some set 
of test vectors. In doing so, we are able to identify suitable coarse-grid points as those from which it is 
possible to build a high-quality LS interpolation to its neighbors. Next, we describe the idea of measuring 
strength between neighbors via algebraic distances in detail and then discuss how we use this measure in a 
CR coarsening algorithm and in computing interpolation. 

4.1. Strength of connection by algebraic distances. The basic idea of the algebraic distances 
approach to measuring strength of connection is to construct a row of least squares interpolation, pi, for 
each fine variable i € f2 and various choices of its interpolatory set, Cf. 

CS(p t ) = » K f v\ K) - ~r\ K) - J2 (Pi)ivf ) ^ min, (4.1) 

and then using the associated values of £S to define weights quantifying the connectivity among variables. 

With a given matrix A, associate a connected graph Q = (V,£), where V and £ are the sets of vertices 
and edges, with n = |V| (cardinality of V). Here, an edge (i, j) G E (A)ij ^ 0. Let Gd(Vd,£d) denote 

the graph of the matrix A d and define Gd,i{Vi,£i) as the subgraph associated with the ith vertex and its 
algebraic neighbors: 

V i :={j\{A d ) ij ^Q} and % := (A% ^ 0}. (4.2) 

In its simplest form, the notion of algebraic distances is straightforward: For a given search depth, d, 
set of test vectors, {v^\ ...,v^}, and fine variable ieSJ compute 

for all j e Vi, (4.3) 




Ek ( (k) l (k) (V 

where pij is the minimizer of the least squares problem to a single variable j. Here, although the measure 
is able to determine the coupling between any given two points, we limit its use to local neighborhoods, 
j € Vi. This simplification coupled with the idea of deriving strength according to an algebraic distance 
measure based on simple one-point (one-sided) interpolation, allows us to control the complexity of the 
algorithm. More generally, the algebraic distance measure can be defined from the LS functional obtained 
by building interpolation from sets of points, thereby providing a local a posteriori measure of the quality 
of the interpolatory set d. 

Remark 4.1. Using the given weighted strength graph in a coloring algorithm to aggregate unknowns 
or in a greedy algorithm to define a matching in the graph results in a so-called plain aggregation scheme, 
which in turn can be used as a solver within an AMLI cycle \BCKZll\ . In BBB+ 1 df , such a scheme 



was 



developed for finding steady state vectors of Markov chains. A similar technique, which chooses aggregates 
using a greedy strategy based on a local stability measure is developed and analyzed in \BCZ11 



4.2. Compatible relaxation coarsening and algebraic distances. In choosing C, we integrate 
the simplified variant of the algebraic distance notion of strength of connection based on one-sided interpo- 
lation (4.3) into the CR-based coarsening algorithm developed in [BF10] . The notion of algebraic distances 



is used to form a subgraph of the graph of the matrix A , d = 1, 2, . . . Specifically, the algebraic distance 



between any two adjacent vertices in the graph, Qd, of A d is computed using (4.3). Then, for each vertex, 
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i € F, we remove edges adjacent to i with small weights relative to the largest weight of all edges connected 
to i: 

V M = F ; £ M = {(i,j)\i,j £ F and r i} ;> 9 ad vaasLr lk }, (4.4) 

k 

with d a d G (0, 1). This in turn produces the graph of strongly connected vertices, A^ d (V»,£x). Note that 
by definition, the strength graph is restricted to vertices i G F so that it can be used in subsequent CR 
coarsening stages. 

The strength graph, A4 d , is then passed to a coloring algorithm, as in the classical AMG approach 
outlined in Section [2j where as mentioned earlier, coarse points are selected based on their number of 
strongly connected neighbors. Recall that, in the classical approach strong couplings are defined in terms of 



Si and Sf as in Equation (2.1 1. In our approach, the sets Si and Sf containing the information on strong 
coupling are instead derived from the algebraic distance based connectivity graph A4 d . 

A rough description of our overall CR coarsening approach is described by Algorithm 1; for additional 



specific details of the CR algorithm we refer the reader to BF10 , in which this scheme was developed for 
diffusion problems similar to those we consider here. 



Algorithm 1 compatible_relaxation {Computes C using Compatible Relaxation} 

Input: Co {Co = is allowed}. 
Output: C 
Initialize C = Cq 
Initialize Z = Q \ C 

Perform v CR iterations with constant u° 
while Pf > 6 do 

Z = {iefl\C :<7i> tol} 

C = C U { independent set of Z } based on jVC' 
Perform v CR iterations with constant u° 
end while 



Remark 4.2. Various choices of the candidate set measure, o~i, used in determining potential C-points 
have been studied in the literature fBraOO, Liv04 , BFlfJj ; our 's follows the definition in JBF1Q 1, namely 



We mention that, in practice, this measure gives the best overall results for smaller values of v, say v = 5. 
Additional discussions and extensive testing of this approach are found in fBFlCtf . 

Remark 4.3. The matrix A d , d > 1, is never formed in practice - it is straightforward to reconstruct 
the connectivity graph of A d using local operations involving only the graph of A. 

4.3. Denning interpolation via algebraic distances. Given a set of coarse variables C and a set of 
r-smoooth test vectors, {v^\ v^}, interpolation of a given caliber, c, is constructed via the least squares 
functional defined in ( |4.1[ ). More specifically, we consider all possible sets of C-points, W, up to a given 
cardinality as prescribed by parameter, c, in the dig-ring coarse point neighborhood of a given F-point, i, 
defined as 



Af dis ,i := Cn V dl 



(4.5) 



where the connectivity in Gd LS ,i is defined as in (4.2). Thus, the sets of possible interpolatory points can be 
written as 

W := {W\ W C Af dLsA and \W\ < c}. 
By sampling this set we find the minimizer of the least squares functional: for each i G F 



Ci = arg min £S w (Pi), 
wew 
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where, 



k 

CS w {Pi) = },u K 

K=l 



,(«) 



(«) 



and denotes the minimizer of the least squares problem (3.2) for the given set W. Thus, we must compute 
Pi and evaluate CSw{Pi) for all possible choices of W G W. The row of interpolation, pi, is then chosen as 
the one that minimizes CSw{Pi)- 

An additional detail of our approach is the penclization of large interpolatory sets. It is easily shown that 
for two sets W C W", their corresponding minimal least squares functional values fulfill CSw' > CSw n - In 
order to keep the interpolation operator as sparse as possible, we require that the least squares functional 
is reduced by a certain factor when increasing the cardinality of W . That is, a new set of points W" is 
preferred over a set of points W with \W"\ > \W'\ if 



CS 



< {cs w ,)^ w "^ w '\) 



Based on numerical experience we usually choose 7 = 1.5 which tends to produce accurate and sparse 
interpolation operators for a large class of problems. 

Remark 4.4. The above sampling of all possible combinations of interpolatory sets with cardinality 
up to some given caliber is one of many possible strategies for selecting Ci. In our experience, such an 
exhaustive search is rarely necessary and in most cases scanning a small number of possibilities using some 
systematic strategy that incorporates the algebraic distance strength measure based on one-sided interpolation 
is sufficient. 

5. Numerical Results. In this section, to demonstrate the effectiveness of the algebraic distances 
measure of strength of connection when combined with CR coarsening and a greedy approach for constructing 
LS interpolation, we present tests of the approach applied to a variety of 2D anisotropic diffusion problems 
discretized on a (N + 1) x (N + 1) uniform grid. 



5.1. Model problem and its discretization. 

two-dimensional diffusion operator 



We consider a finite-difference discretization of the 



Cu = V • /CVu 



with anisotropic diffusion coefficient 



( cos a 
\ sin a 



— sm a 
cos a 



1 

e 



cos a 
- sin a 



sm a 
cos a 



where < e < 1 and < a < 2tt. Changing variables 



cos a sm a 
- sin a cos a 



yields strong connections aligned with direction £: 

£u(£,,T]) = Utf + eu vv . 
Its equivalent formulation in (x, y) coordinates is given by 

Cu(x, y) = au xx + bu xy + cu yy , 
where a = cos 2 a + e sin 2 a, b = (1 — e) sin 2a, and c = sin 2 a + e cos 2 a. 



(5.1) 



(5.2) 



(5.3) 



(5.4) 



(5.5) 



In formulating a finite-difference discretization of (5.5 1, we consider a standard five-point discretization 
for the Laplacian term and then define the discretization of the u xy term using intrinsic strength of con- 
nections. In the a = 7r/4 case, for example, this amounts to using lower-left and upper-right neighbors and 
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avoiding lower-right and upper- left ones in denning the discretization for each fine-grid point i G fl. The 
stencil for u xy and this choice of a is then given by 

1 /o -1 1 \ 1 / o 0\ 

(5.6) 

The overall stencil thus includes seven nonzero values. In our numerical experiments, we consider also the 
worst case scenario, that is, the case where the discretization for u xy is defined along the weakest connection, 
r\. For instance, the appropriate stencil for a = n/A is a poor choice for a = — 7r/4 as directions £ and ry, 




as in (5.3 1, interchange for these choices of a. Further, for this value of a and taking e = .1, the resulting 




system matrix has stencil 

S A = 

and, thus, is not an M-matrix. Here, some of the off-diagonal entries in A are positive and, hence, the 
heuristics motivating the classical definition of strength of connection are not applicable. We consider this 
extreme case, although it is unlikely to arise in practice, to demonstrate the robustness of our approach 
in choosing the right coarsening strategy for the targeted anisotropic problems. We mention that on a 
structured grid, our chosen seven-point discretization is equivalent to the finite element discretization of the 
same elliptic boundary value problem. 

5.2. Formulating a two-level coarsening algorithm. Our aim is to study the robustness of the 
algebraic distance notion of strength of connection for grid and non-grid aligned anisotropy. We focus our 
tests on the seven-point finite-difference discretization introduced earlier in this section for various values of 
anisotropy angle, a, anisotropy coefficient, e. 



In all tests, coarse grids are chosen using the compatible relaxation approach discussed in Section 4.2 and 
then interpolation is computed using the greedy LS approach given in Section |4.3| The stopping tolerance 
for CR steps is set as 6 = 0.7 and the number of CR sweeps is set as v = 5, so that the algorithm terminates 



when pf in (3.4) is below this value, and the tolerance for the candidate set measure is taken as tol = 1 — pt. 
The larger choice of stopping tolerance allows more aggressive coarsening whereas our choice of tol ensures 
that points are added to the candidate set more sparingly at later stages of the CR coarsening scheme. We fix 
the strength of connection parameter used in forming the strength graph at 9 a d = -5 and vary the the graph 



distance, d, used to define the graph Qd, from which Aid in (4.4 1 is constructed. We note that generally the 



overall quality of the grids the algorithm produces depends only mildly on the choice of 9 a d- 

The search depth, used in defining the greedy approach for choosing LS interpolation as in (4.5), is 
fixed as cIls — d + 2. Here, taking a larger value of the search depth than the coarsening depth allows the 
approach to scan a larger number of possible interpolatory sets in forming long-range interpolation (whenever 
the problem requires it). In this way, the LS scheme for constructing interpolation is able to better follow 
a wider range of values of the anisotropy angle, a. Eight test vectors are used to construct the modified 



LS interpolation given in (3.2), seven are computed by applying 40 iterations of Gauss Seidel to distinct 
random initial guesses and the eighth is the constant vector also relaxed with 40 iterations of Gauss Seidel. 
We mention that the number of test vectors and especially the amount of relaxation applied to them can be 



sufficiently reduced by replacing a one-grid procedure by a multilevel bootstrap cycling scheme BBKLIf 
However, the main focus of this paper is the study of the performance of the algebraic distances measure of 
strength of connections for the targeted anisotropic problems; the development of a multilevel algorithm is 
subject of on-going research. 

When presenting results of the two-grid solver constructed by our algorithm we use two pre- and post- 
Gauss Seidel relaxation sweeps on the fine grid and a direct solver on the coarse grid. Here, we use two pre- 
and post- smoothing steps because our algorithm generally produces aggressive coarsening. The estimates 
of the asymptotic convergence rates are computed as follows: 

|e"|U 



9 lie"- 1 ! 
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where e n is the error after r\ = 100 MG iterations, applied to the homogenous system starting with random 
initial approximations. We also report the operator complexity ratio 

nnz(A) + nnz(A c ) 
nnz(A) 

\C\ 

and the coarsening factor j g — — . 

5.3. Choosing coarse variables and interpolation. In this section, we present the choices of coarse 
grids our algorithm makes when applied to the anisotropic problem for e = 10~ 10 and various choices of 



the anisotropic angle a in (5.2). In the plots in Figures 5.1 - 5.3 the black lines depict the interpolation 
pattern for each _F-point (denoted by the smaller circles) from its neighboring C-points (denoted by the 
larger circles). We note that for the d = 1 tests in Figure |5~Tj the coarsening and interpolation pattern follow 
the anisotropy when the choice of the coarse grids allow it, i.e., for a — and a = n/A. This is not the 
case for the a — —ir/A and a = ir/8 cases, for which the direction of strongest coupling is not captured by 
the connectivity of the fine-level discretization of the problem. Hence, for these problems it is not possible 
to form a strength matrix from immediate algebraic neighbors which produces a coarse grid that allows the 
interpolation pattern to follow the anisotropy. These observations in fact led to the idea of using the graph 
Gd of A d , d > 1, to form the strength matrix. Overall, we see that even for these cases, the algebraic distance 
strength measure leads to least squares interpolation that follows the general direction of anisotropy as much 
as possible when coarsening is done using Q\ to form the strength graph and set of coarse grid points. 



Another interesting observation here is that for the d — 2 results in Figure 5.2 by using the graph of 
A 2 to form the algebraic-distance-based strength graph, the coloring algorithm is now able to coarsen in 
the direction of anisotropy for the a — — 7r/4 case. This choice also allows a more aggressive coarsening for 
other values of a while maintaining, for a — and a = 7r/4, the characteristic directions induced by the 
anisotropy. 

A particularly difficult choice of a occurs when the anisotropy direction is nearly aligned with a grid 
direction, for example, taking a ~ for our diffusion problems on a regular 2D grid. For these cases, a 
longer-range interpolation and an extended search depth for coarse-grid candidates is needed to accurately 
capture the anisotropy (more so for directions closer to the grid lines). To demonstrate this phenomenon, 



we include the results for a = tt/8 in Figure 5.3 Here, we see that when using a larger value of d in (4.4) to 



define Qd and also a larger value of the search depth (d^s in (4.5)), the LS procedure is able to generate an 



interpolation operator that accurately follows the direction of the anisotropy. 

5.4. The coarse-level operator. One of the interesting deliverables of the algorithm, in particular of 
its implementation of the compatible relaxation and the algebraic distances, is the pattern of the resulting 



coarse-grid stencil. Discretization involving (5.6) favors a = tt/A and with the same argument makes the 
worst possible discretization for a = — it /A for which employment of upper- left and lower- right grid-point 
neighbors in discretization of d x y wo uld be beneficial. We now consider both cases a = ±7r/4 and the seven 



point discretization, employing ( |5.6[ ). For both cases, we assume e = 10 10 , d = 2, and dLS = 4. 

We first confirm that for a = 7r/4, the coarse-grid operator A c = P t AP preserves the intrinsic strength of 



connections inherited from the fine-grid operator A. A typical example of its stencil is given in Fig 5.4(a) (the 
details of configurations depend on the coarsening pattern in the neighborhood of the considered coarse-grid 
equation) . 



The results for the more challenging a = —tt/A case are provided in Figure 5.4(b) Here, we observe that 
although the discretization on the fine grid does not follow the anisotropy whatsoever, the non-zero pattern 
of the coarse-grid operator correctly aligns with the direction of anisotropy. This result demonstrates the 
ability of the algorithm to overcome, if needed, the disadvantage of a poorly chosen fine-grid discretization 
and regain a more favorable discretization on the first coarse grid. Further, the results for the a — tt/A 
indicate, that all consecutive coarse grids (though not employed in our two-level algorithm) are likely to 
maintain a similar favorable discretization, that too accurately reflects the anisotropy. 



Coefficients of the coarse-grid stencils, presented in Fig 5.4 are given next. Here corresponds to 
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Fig. 5.1. Coarse grids and interpolation patterns for h = 1/32 for various choices of a, using the graph of A, i.e., d = 1 
and d^s = 3, to define the strength matrix. 



a = ir/4 (entries denoted by * are negligible, with absolute values below 10 n ) 



* -0.166 \ 



9+ - 

D cg — 



0.332 



-0.166 
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5 10 15 20 25 30 

(a) a = tt/8, p = .29, p/ = .59, -y a = 1.805, -y g = .295 




5 10 15 20 25 30 

(b) a = tt/8, p = .40, p f = .76, 7 = 1.543, -y g = .257 



Fig. 5.3. Coarse grids and interpolation patterns for h = 1/32 with a = n/8 using the graph of A 2 , i.e., d = 2 and 
di,s = 4 (left) and A 4 , i.e., d = 4 and d^g = 6 (right). 



30- 



25- 



20- 



15- 



10 



5-' 




30- 



25 • 



20- 



15-' 



10 



5- ■ 




5 10 15 20 25 30 

(a) Coarse-grid equation pattern for a = n/4. 



5 10 15 20 25 30 

(b) Coarse-grid equation pattern for a = — 7r/4 



Fig. 5.4. Non-zero pattern of one of the stencils of the coarse-grid operator A c , centered at i £ C and connected with 
j S C such that (A c )ij 0. The smaller dots in the graph are all other coarse-grid variables. 



In both stencils, all entries are rounded to the nearest hundredth. 

The distribution of the stencils' coefficients further illustrates the ability of the algebraic distance based strength 
of connections measure to choose the correct coarse-grid points for problems with anisotropic coefficients; in both 
cases, the non-zero sparsity pattern and dominant coefficients of the resulting coarse-grid operators follow the direction 
of anisotropy. 

Remark 5.1. We recall that a larger graph distance d — 2 and search depth dLS = 4 were required to obtain our 
promising results for the a = n/4 case and note the additional fill-in of the resulting coarse-level operator in this case, 
as seen in Figure \h\4\ b. Generally, as the direction of the non-grid aligned anisotropy aligns itself more closely with 
the grid, the values of d and, hence, dLS, must be increased for the approach to appropriately capture the anisotropy. 
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This, in turn, increases the fill-in of the coarse-level operator, making it difficult to recursively coarsen the equations 
in a systematic way and maintain low grid and operator complexities. We mention that any method for constructing 
the long-range interpolation required by such problems will have to be designed with this issue in mind. 

5.5. Two-level convergence. We conclude our experiments with tests of the proposed AMG setup algo- 



rithm applied to (5.1 1 for various choices of the anisotropy angle, a, the anisotropy coefficient, e, and the mesh 
spacing h = 1/N. In the following tables, the asymptotic convergence rates, p, of the two-grid solver produced by 
our setup algorithm are reported, along with the corresponding coarsening factors 7 9 and operator complexity ratios 
7 . Here, we observe a slight dependence of the computed convergence rates and grid and operator complexities on 
the problem parameters, e, a, and h. This dependence on h is restricted mostly to the non-grid aligned cases, with 
the exception a = and e — .1, where we see a slight increase in p as the problem size is increased from N = 64 to 
N = 128. Moreover, in all cases, the convergence rates and complexities are uniformly bounded with respect to e and 
a for fixed h. We note in addition that these results are promising when considering that all tests were performed 
with the same strength parameter 6 a d = -5. In fact, all parameters in the setup algorithm were fixed, suggesting that 
the individual components of the setup are robust for the targeted anisotropic problems, even those leading to non 
M-matrix systems as in the a — —n/4,n/8 cases. Further, we note that the setup handles the isotropic case when 
a = and e = 1 with similar efficiency, producing a two-grid method with convergence rate p — .28 and complexities 
7 9 = 25 and -y = 1.6 for N = 32, 64, 128. 



a 


6 = .1 


tt/4 


.10 (.35,1.5) / .22 (.30,1.5) / .22 (.27,1.5) 


-tt/4 


.31 (.35,1.7) / .36 (.34,1.7) / .48 (.32,1.7) 


tt/8 


.32 (.27,1.4) / .39 (.24,1.5) / .35 (.25,1.5) 





.19 (.34,1.6) / .20 (.37,1.7) / .24 (.38,1.7) 




a 


e = .0001 


vr/4 


.26 (.35,1.5) / .26 (.34,1.5) / .23 (.36,1.5) 


-tt/4 


.28 (.46,1.9) / .33 (.45,1.9) / .38 (.41,1.9) 


tt/8 


.30 (.43,1.8) / .48 (.38,1.8) / .51 (.36,1.8) 





.05 (.34,1.6) / .06 (.35,1.6) / .06 (.38,1.7) 




a 


e = 


tt/4 


.06 (.33,1.4) / .06 (.34,1.4) / .06 (.36,1.5) 


-tt/4 


.28 (.46,1.9) / .35 (.43,1.9) / .37 (.41,1.9) 


tt/8 


.30 (.43,1.8) / .49 (.38,1.8) / .52 (.36,1.8) 





.05 (.34,1.3) / .06 (.35,1.3) / .06 (.38,1.4) 



Table 5.1 

Approximate asymptotic convergence rates for the seven point FD Anisotropic Laplace problem with Dirichlet boundary 
conditions for various choices of a, e and h. Here, our proposed setup is applied with d = 2 and d^g = 4. The reported results 
correspond to: p (7 fl ,7o) for n = 32 2 / n = 64 2 / n = 128 2 . 



6. Concluding remarks. The LS functional gives a flexible and robust tool for measuring AMG strength of 
connectivity via algebraic distances: between pairs of points to define a strength graph used to choose coarse points; 
and among sets of points for determining interpolatory sets. The proposed coarsening approach combining algebraic 
distances, compatible relaxation, and least squares interpolation provides an effective scheme for the non-grid aligned 
anisotropic diffusion problems considered. The approach chooses suitable coarse-grid variables and prolongation 
operators for a wide range of anisotropies, without the need for parameter tuning. Moreover, even when the initial 
fine-grid discretization is chosen in the direction opposite to the one defined by the anisotropy (as in the a — — 7r/4 
case), the method constructs a suitable interpolation operator and, further, produces a coarse-grid operator which 
better captures the anistropy. While not the focus of this paper, the promising results obtained by our approach 
suggest that its extension from two to many levels should be effective. Indeed, it is natural to assume that if the 
first coarse grid system produced by our scheme is consistent with a suitable discretization, regardless of the initial 
discretization, then the multilevel scheme constructed from it will yield an effective solver. As noted earlier, the 
main challenge faced in extending our approach to a multilevel one for the targeted anisotropic problems is that of 



Algebraic distances in AMG 



15 



designing an algorithm capable of constructing long-range interpolation as needed to capture general anisotropies 
and at the same time maintains low grid and operator complexities. 
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